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Until recently it has been impossible to 
accurately determine the roots of poly- 
nomials of high degree, even for polyno- 
mials derived from the Z transform of 
time series where the dynamic range of 
the coefficients is generally less than 100 
dB. In a companion paper, two new pro- 
grams for solving such polynomials were 
discussed and applied to signature anal- 
ysis of one-sided time series [1], We 
present here another technique, that of 
root projection (RP), together with a 
Gram-Schmidt method for implementing 
it on vectors of large dimension. This 
technique utilizes the roots of the Z 
transform of a one-sided time series to 
construct a weighted least squares modi- 



fication of the time series whose Z 
transform has an appropriately modified 
root distribution. Such a modification 
can be employed in a manner which is 
very useful for filtering and deconvolu- 
tion applications [2]. Examples given 
here include the use of boundary root 
projection for front end noise reduction 
and a generalization of Prony's method. 
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1. Introduction 



Let y be a complex number. We define 'y' to be 
the "geometric sequence" vector with components 
"yi^^y", /j=0,..., N. Then given a time series of 
length AT-i-l represented by the vector a:N, 
a(y) = '^^a ; and the condition that y be a root of 
a(y) is: 



y'a =0. 



(1) 



Given a (possibly complex) time series a:N and any 
set of complex roots of eq (1), y\,,„, yu, M^N,-we 
shall construct a linear projection to modify a:N to 
another time series b:N so that b(ya) = 0, a = l,..., 
M, and so that 



ia-bYGia-b)*, 



(2) 



where "*" means complex conjugate, is a minimum 
for any positive definite Hermitian weighting ma- 
trix G. In other words b:N is the series closest to 
a:N in a weighted least squares sense, which has a 
chosen set of complex numbers among the roots of 
its Y transform. We refer to this process as root 
projection (RP). 

In a slightly more general context consider a 
(possibly complex) time series a:N, a set of com- 
plex numbers y\,,.., yu with M^A^, and a set of 
complex values vi,,.,, vm. We construct a time series 
b:N whose difference from a:N is minimal among 
those series with b (ya) -Va, a = 1,..., M. 

Because of the slight complication due to the 
complex numbers involved, we shall set and solve 
the above minimization problem using Lagrangian 
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multipliers. We want to find the time series b:N 
such that: 

M M -I 

2:A*(i:^„ -v.)+^K{x*-^yt -wj) .(3) 

a=l a=l J 

Using index notation and differentiating eq (3) 
with respect ;c/* and A*, we have: 






(4) 



S&nmg y^={y,f, yi, = {y*f, and G,7' =(G "'),;; 
solving for b, in the first of eqs (4), substituting into 
the second equation and solving for A^: 



(5) 



Ap = 2 pj {ujy^ - Va), where 

QaP =yakGu^1i 

whence: 

bi =ai -ajy^Q^y%,G;;,l + v,jQ^y%,G;;u^ . (6) 

To make the connection to linear vector space 
theory, we start with complex N space, (£", having 
an inner product given by the positive-definite Her- 
mitian matrix G: 



(a,b)=a^G''b*, 



(7) 



and complex M space, (£^, having an inner product 
given by the Hermitian matrix Q "' (assumed, here, 
to be non-singular) with an inner product defini- 
tion similar to eq (7) and with vectors such as v 
(with components u„, a = 1,..., M) [3]. We define 
the mapping Y from (£" to ©^ by: 

Y=[Yu =y^; a = 1,..., M, i = 1,..., N]. (8) 

Y has an adjoint mapping Y^, from (S," to S" given 
by the condition (Y^vyG^a* = v\Q-Y(Ya)*, 
from which 

Y^ = G-^Y*'^Q-K (9) 

In matrix notation eq (6) becomes: 

b = {I-P)a+Y^v 

(10) 
P = Y^Y,P^=P, 



the projection property, P^=P, following directly 
from eq (6). It also follows directly from eq (6) that 
the range of /— f is contained in the null space of 
Y^, while the range of P is contained in the range 
of Y^. Since both I—P and P as well as the null 
space of y and the range of Y^ decompose (£" (the 
latter being an orthogonal decomposition), I—P 
and P are orthogonal mappings onto the respective 
spaces. The mapping / —P is, then, the desired root 
projection.' 



2. A Gram-Schmidt Algorithm for Root 
Projection 

In order to carry out root projection, we need an 
orthogonal basis for the range of Y^ in the unitary 
space with metric G . From eq (8) we see that this 
space is spanned by vectors of the form G~^y^. For 



any two such vectors: 



t-^ 



(G-'y:rG'-(G-'Tpr=(G-^„r(G~^,r, 

(11) 



where 

(G-if=G-K 



(12) 



Thus, an orthonormal basis chosen from the vec- 
tors G ""y^ by applying the Gram-Schmidt process 
using the ordinary unitary inner product will deter- 
mine an orthogonal basis for the range of Y^ in the 
unitary space with metric G by multiplying each 
vector by G "2. However, when only a limited num- 
ber of vectors need to be projected, the projection 
can be carried out more efficiently by: 

i) transforming the vector a to G^a 
ii) projecting Gia in the Euclidean norm us- 
ing the modified basis G "y^ 
Hi) transforming the result back by multiplica- 
tion with G "2 . (13) 

The root projection process is especially simple 
when the set of ya's is closed under complex conju- 
gation (i.e., if ya belongs to the set, so does y*) for 
the case of simple time or frequency weighting. Be- 
cause of the closure under complex conjugation, 
the y^ vectors can be replaced by the vectors con- 
sisting of their real and imaginary parts, yi(y^) and 



' A concise exposition of the linear algebra required for root 
projection, with G =1 can also be given in terms of the QR 
decomposition [4] (although we use the "transposed" form of 
that in [4]). 
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S(5v ). By the time weighting case we mean that 
G~^ and (3~2 are diagonal (and, therefore, positive 
and real), so that no unitary transformation is re- 
quired. In that case each of the rows is multiplied 
by the same factor. Thus, the basis for the range of 
Y^ can be obtained by applying the real form of the 
Gram-Schmidt process to the vectors dt(G~2y^) 
and ^(G ~y7). The root projection of a real time 
series, outlined in eq (13), will, in this case, also be 
real. 

In the case of frequency weighting we mean 
scalar weights applied to the DFT components of a 
time series. The DFT process is, itself, a unitary 
transformation, and the orthogonalization proces^ 
can be made real by the trick of putting 91(5* )v'2 
in place_of the DFT component a* and putting 
^(at)\/2 in place of the DFT component oT*. The 
real dot product of these series is the same as the 
unitary dot product of the DFT's. A diagonalized 
frequency weighting metric can then be applied to 
these transformed DFT series and the diagonaliza- 
tion process carried out. If projection is to be car- 
ried out on a real time series, it can be done in 
these transformed coordinates, otherwise, the new 
orthogonal DFT coordinates have to be recon- 
structed before doing the projection. Of course, the 
projected DFT series has to be reconverted to time 
series form after projection. 

The real form of root projection can also be ob- 
tained by using QR algorithms such as given in 
LINPACK [4]. However, the implementation of 
the QR algorithms in LINPACK is CPU core con- 
suming for large N and does not easily allow stor- 
age of the orthogonal matrix. These disadvantages 
are overcome using the Gram-Schmidt method 
given above. Employing a form of the Gram- 
Schmidt algorithm suggested by G.W. Stewart 
[5,6], a dynamically dimensioned Fortran 77 code 
was constructed allowing efficient use of a user se- 
lected buffer with the option of saving the orthogo- 
nalized basis vectors on disc for rapid repeated 
projection [7]. For ya with modulus greater than 
one, projection was carried out using the root re- 
ciprocal to avoid overflow. In this case the series 
starts &\.ya~^ and goes up to 1, a procedure equiva- 
lent to reversing the series to be projected. 



3. Examples 

We have principally used root projection for de- 
convolution. That topic is discussed in the third pa- 
per in this series [2]. Here we give four illustrative 
examples of root projection. The fourth example. 



applying root projection to Prony's method may 
have some practical use. However, a detailed study 
has not been conducted. 

1. Figure 1 shows a normalized 101 point Gaussian 
with 60 dB dynamic range (ratio of center to edge 
points). In figure 2 we show the 800 point result of 
convolving the first 700 points of the experimental 
waveform of figure 2.4a in [1] with the Gaussian of 
figure 1. The noise pattern shown in figure 3 re- 
sults from rounding off the time series in figure 2 
to 8 bit accuracy. The rounded-off time series is 
not shown. The standard deviation of this noise is 
6.02x10"". 

A basis for all 800 point time series can be 
formed from the waveform of figure 2 together 
with the 799 geometric root vectors formed from 
the roots of the Y transform of that waveform. The 
noise vector of figure 3 must then be a linear com- 
bination of all these vectors, while the true signal 
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Figure 1. Normalized 101 point Gaussian with 60 dB dynamic 
range. 
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Figui^ 2. Convolution of the first 700 points of the time series 
in figure 2.4a of reference [1] with the 101 point time series of 
figure 1. 



vector of figure 2, which is one of the basis ele- 
ments, is orthogonal to all of the 799 geometric root 
vectors. If we use the geometric root vectors built 
from the Gaussian of figure 1 and apply root projec- 
tion to the rounded-off series built from the series 
in figure 3, the noise of the resultant series will be 
reduced in magnitude, since it is orthogonal to the 
100 geometric root vectors formed from the roots of 
the y transform of the Gaussian. The standard devi- 
ation of that noise is 5.6 x lO""*, almost exactly the 
theoretical estimate of (5)2 expected for Gaussian 
noise. If the 799 geometric root vectors from both 
the Y transforms of the Gaussian and the waveform 
of figure 2.4a in [1] are used for the projection, then 
the error appears as in figure 4— a mini-image of 
the correct time series— with a standard deviation 
of 2.0x10"^ (again approximately (800) "2 of the 
original noise deviation). This time, of course, the 
noise is biased with a mean value of 4.9 x 10"'. This 
curve also indicates the precision of the root finding 
and Gram-Schmidt routines used in these calcula- 
tions. 
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Figure 3. 8-bit roundoff noise for the series in figure 2. 
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Figure 4. Error from projecting the series in figure 2 plus noise 
through all roots of the Y transform of the series in figure 2. 
Note that the remaining error is proportional to the original 



2. The use of an extensive causal boundary root set 
for front end filtering is shown in figures 5a and 5b. 
To the Gaussian filtered curve of figure 2,3 a in [1] 
we add the - 40 dB noise distribution shown in fig- 
ure 5a and apply root filtermg using only the causal 
boundary roots from figure 2.1b of [1]. The noise 
after filtering showing extensive front end reductiott 
is shown in figure 5b. 

3. A more general example of root projection em- 
ploying both time and frequency weighting is af- 
forded by constructing an approximation to an 
optimal maximal-ripple lowpass filter. Figures 6 
shows such a filter of 151 elements designed using 
the Remez exchange algorithm of McClellan et al. 
[8]. In this case the passband was set for 20% of the 
Nyquist frequency and the stopband for frequencies 
over 33% of the Nyquist frequency. 

We start with a delta Series shifted So that the 
value 1,0 lies at position #76 in the Center of the 
time interval. To use unweighted projection and 101 
roots evenly spaced from - 60 ' to +60 ° on the unit 
circle would produce the familiar series obtained 
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Figure 5a. Noise added to series in figure 2.3a of reference [1]. 
337 



Volume 96, Number 3, May- June 1991 

Journal of Research of the National Institute of Standards and Technology 



U. lb- 


— , — , — , — ■ — !_, , — , — , — t — 1 — , — . — , — !_, , , — , — t — , — ,_ 

1 1 


_, — , — !_, — , — , — , — I — , — , — , — , 


1, 




■ 


0.09- 


1 


, 1 
1 1 , 






■ 


^ 0.03- 




I 1 1 

1 1 


1 ll 




1 


=D 




1 




lllllilllllllllLIf 


h- 






jl l'{lllillllllllltllll 


1 — 1 
Ql 

^ -0.03- 


^ 




i 1 

j 

1 1 1 


1 ■ 


-0.09- 




i 




1 
1 


j 


-n 1 q- 




' 


1 







G.OG 100.00 200.00 300.00 
POINT NUMBER 



400.00 



Figure Sb. Noise after filtering using causal boundary roots of the series in figure 2.3a of 
reference [1] showing extensive front end noise reduction. 



from windowing the DFT of the shifted delta se- 
ries. However, by using time weighted projection 
with the center symmetric time weight shown in fig- 
ure 7 and adding two extra roots near each edge 
of the stopband region, one can produce a filter 
with 120 dB attenuation in the stopband and 
—2.5x10"' leakage in the time domain, but with 
0.11 ripple in the passband. Without time weight- 
ing the filter produced has 64 dB attenuation in the 
stop band and 0.08 ripple in the passband, but with 
— l.Ox 10"^ leakage in the time domain. 

To remove the passband ripple from the 
weighted filter we subtract from it the shifted delta 
series. The resulting series should be close to zero 
in the passband region. We place 31 roots evenly in 
the passband and again place two extra pairs of 
roots near the edge of the band. Applying fre- 
quency weighted projection with a simple square 
frequency weight to hold the stopband attenuation 
in place, we obtain (after re-adding in the shifted 
delta series) the doubly root projected filter shown 
in figure 8, where it is compared with the un- 
weighted doubly projected and optimal filters. The 



weighted series closely resembles the optimal series 
of figure 6 with leakage at the ends of the time 
interval of -8.0x10"^ as opposed to ~2.0xl0~^ 
for the unweighted filter. The leakage for the opti- 
mal filter is ~ 6.0x10"*. The attenuation in the 
stopband is 70 dB and the ripple in the passband is 
1.0x10"'' for the weighted filter while the 
unweighted filter has 36 dB attenuation and 
8.0 X 10"" ripple. The optimal filter has 116 dB at- 
tenuation and 5.0 x 10"^ ripple. The spectra for the 
three filters are compared in figure 9. 

The root patterns in the stopband for the opti- 
mal filter are shown in figure 10. As can be seen, 
they also cluster near the edge of the stopband, but 
they are not evenly spaced throughout the stop- 
band. Similarly, the roots of the optimal filter mi- 
nus the shifted delta series, shown in figure 11, 
exhibit the same uneven spacing and clustering in 
the passband region. 

4. Root projection can be used to provide a 
"global" least squares generalization to the Prony 
method. Prony's method is applied to the case of 
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Figure 6. Optimal maximal ripple filter of 151 elements. 

N + 1 data points which one wants to represent as 
the sum of n decaying exponentials (Prony takes 
N + 1 even and n = (N + 1)/2): 



F, = ^jyf,k=0,...,N. 



(14) 



Here yj is thought of as a complex number of mod- 
ulus less than 1, yj = e^^, real()^) < 0, and T a sam- 
pling interval (See, e.g., [9]). Eq (14) can be 
rewritten, taking advantage of the cyclotomic poly- 
nomial expression, as 

F(y)=lF,y^=lA ^\^jy\ ■ (15) 

k=o ;=i J- yjy 

Combining the terms in eq (15) together to form 
a single fraction gives 
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Figure 7. Time weights for approximating an optimal filter. 

F(y) = [a^(1 -y2yyil -yny){l - (yiyf*') + ... 

+An(l-y,yyil-yn-iy)[l-(ynyr')]/ 

[(l-yiy)-(l-yny)]. (16) 

Equation 16, then, has the form 



^^ Q(y) 



(17) 



where Q(y) is a polynomial of degree n whose 
roots are (yj)'^, and P(y)+y'^*^R(y) is a polyno- 
mial of degree N+n whose N + l-n coefficients 
from that of y" to that of y" are zero. The roots of 
Q(y) are referred to as the poles of F. Conversely 
one can show that if eq (17) holds with the corre- 
sponding coefficients zero, then setting 
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Figure 8. Comparative time plot showing the effects of using weighted projection for approximating an optimal filter. 
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Figure 9. Comparative frequency plot showing the effects of using weighted projection for approximat- 
ing an optimal filter. 
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Figure 10. Roots of the Y transform of the optimal filter shown in figure 
6 showing rcxits in the stopband. 
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Figure 11. Roots of the Y transform of the optimal filter of figure 6 
minus the shifted delta series showing roots in the pass band. 
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e(y)=n(i-M) 



y=i 






(18) 



yields the expression 14. 
The goal, then, is to solve the algebraic equation: 



e(yF(y)=/'(y)+r^'i?(y), 



(19) 



where P, Q, and R are all unknown. Only the in- 
teger n is given. This is generally a very ill-condi- 
tioned process, but when thej'y's are real or occur 
in conjugate pairs, the following root projection 
procedure is often successful: 

a. Use an initial estimate for the coefficients, c,-, of 
P +y^^^R by setting d = ki for i =0,..., n - 1, c, =0 
for i=n,..., N, and c = K2 for i=N + l,..., N+n; 
and use for the diagonal weighting matrix, W, the 
series wi with h'; = Ai for 1 = 1,..., n, wi = l for 
l=n + l,..., N + 1, and H'; = A2 for l=N + l,..., 
N+n + 1, where ki, K2, Ai, and A2 are user selected 
values. 

b. Project the series made from P +y'^'^^R into the 
range of F using W and the geometric row vectors 
made from the roots of F(y). The projected poly- 
nomial will have negligible values for the coeffi- 
cients between n + 1 and N + 1 (depending on the 
noise in the data) and will be divisible by F(y). 
Further, P(y) can be read off from the projection 
and Q(y) found by simple division (using FFT 
methods, for instance, as discussed in [2]). 

As long as: i) the projection is stable, ii) the 
projected vector is not zero, and Hi) the values of 
the weights are large enough to produce an answer 
within the noise for similar cases with no data 
noise, then the values of ki, K2, Ai, and A2 have very 
little effect on the calculated j'y's and^i/s. In the 
ordinary Prony case there are as many degrees of 
freedom added in R as are restricted between P 
and R (one degree is an irrelevant multiplicative 
constant since F is expressed as a ratio). 

This method has been tried for several examples 
with both real and complex poles with maximum to 
minimum root amplitude ratios up to 60:1 and with 
minimum root separation down to 0.01. Using dou- 
ble precision data (with Kl = 10^ K2=10'', Ai = 10^, 
and A2 = 10'^) it works well up to about eight poles 
after which the positions of the larger poles start to 
degenerate (due to the ill-conditioning of the prob- 



lem). The representation of the time series coeffi- 
cients of F remains accurate to a relative error of 
about 10"'^ which represents approximately the 
cumulative accuracy of the algorithms involved. Us- 
ing data given to 16 bit precision (around 5 place 
accuracy), the method is able to resolve about 3 
poles (where K■l = 10^ K2 = 10^ Al = 10'^ and 
A2 = 10*) with the positions of the larger poles again 
beginning to degenerate first. 



4. Summary 

The idea of root projection (RP) was introduced 
to permit least-squares modification of a one-sided 
time series allowing a set of given complex num- 
bers to be roots of the Y transform of the time 
series. The general framework of the method was 
presented, and techniques were given to adapt the 
method to the Gram-Schmidt algorithm. For time 
and frequency least-squares weighting, a dynami- 
cally dimensioned form of the Gram-Schmidt al- 
gorithm has been developed to carry out root 
projection. The code has been employed for root 
projection on time series with up to 1600 points 
and achieved better than 12 place accuracy. Illus- 
trative examples of root projection were shown for 
noise reduction and filter construction as well as a 
least-squares extension of Prony's method. 
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